# This file assesses the relationship between perceptions of objective DAs and outcomes conditional on matched pair.
library(dplyr)
library(ggplot2)
library(lme4)
library(lmerTest)
library(fixest)
library(here)
library(xtable)
library(robustbase)
library(estimatr)
library(brms)

load(here("Design", "matches_DA_new.rda"), verbose = TRUE)
load(here("Data", "wrkdat_DA0_new.rda"), verbose = TRUE)
source(here("Design", "nonbimatchingfunctions.R"))

wdat0 <- inner_join(wrkdat_DA0_new %>% select(-geometry), matches_DA_new)

## Add information about which unit is ranked higher versus lower within pair
wdat0 <- wdat0 %>%
  group_by(pair) %>%
  mutate(
    perc_more = rank(vm) - 1,
    social_capital01_rank = rank(social.capital) - 1,
    community_resp01_rank = rank(community.resp01) - 1
  ) %>%
  ungroup()

table(wdat0$perc_more, exclude = c())
table(wdat0$social_capital01_rank, exclude = c())
table(wdat0$community_resp01_rank, exclude = c())

lmer_social_cohesion <- lmer(social.capital01 ~ vm.norm2 + da_prop_vm_20pct_06 + (1 | pair) + (1 | dauid),
  data = wdat0
)
summary(lmer_social_cohesion, corr = FALSE)

lmer_community_efficacy <- lmer(community.resp01 ~ vm.norm2 + da_prop_vm_20pct_06 + (1 | pair) + (1 | dauid), data = wdat0)
summary(lmer_community_efficacy)

social_cohesion_mlm_res <- c(summary(lmer_social_cohesion)$coef[
  "vm.norm2",
  c("Estimate", "Std. Error")
], confint(lmer_social_cohesion, parm = "vm.norm2"))

community_efficacy_mlm_res <- c(summary(lmer_community_efficacy)$coef[
  "vm.norm2",
  c("Estimate", "Std. Error")
], confint(lmer_community_efficacy, parm = "vm.norm2"))

res <- rbind(
  "Social Cohesion" = social_cohesion_mlm_res,
  "Collective Efficacy" = community_efficacy_mlm_res
)
colnames(res)[3:4] <- c("ci1", "ci2")

save(res, file = here("Analysis", "analysis_DA_new.rda"))
